##@file molecules/energy.tcl
# The commands to manipulate the segments
# \verbatim
##
oo::define Charmm method addseg { segid seq } {
    if {[set path [my search $seq]] != ""} {
        #read pdb
        #puts $path
        my sendCmd "bomblev -3"
        my sendCmd "open read unit 10 card name $path
                    read sequ pdb unit 10
                    gene $segid setup
                    close unit 10
                    open read unit 10 card name $path
                    read coor PDB unit 10
                    close unit 10"
        my sendCmd "define missing select hydrogen .or. .not. initialized end"
        my sendCmd "ic fill\nic param\nic build\nhbuild sele hydrogen end"
        my sendCmd "cons fix sele .not.  hydrogen end"
        my do {
            minimize sd 300
            minimize abnr 300
        } until {$ener(ENER) < 0.0 && $ener(GMRS) < 10.0 }
        my do {
            puts "Fix hydrogen"
            minimize sd 300
            minimize abnr 300
            sendCmd "dynamics start verlet nstep 2000 timestp 0.001 nprint 100 -
         iasvel 1 firstt 200 finalt 0 -
         iseed 309234354 -
         inbfrq -1 imgfrq -1 ihbfrq 0 ilbfrq 0 iunwrite -1"
            minimize sd 300
            minimize abnr 300
        } until {$ener(GMRS) < 1 }

        my sendCmd "cons fix sele none end"
        my sendCmd "cons harm force 100 sele .not. hydrogen end"
        my do {
            minimize sd 300
            minimize abnr 300
            sendCmd "dynamics start verlet nstep 2000 timestp 0.001 nprint 500 -
         iasvel 1 firstt 298 finalt 0 -
         iseed 3094354 -
         inbfrq -1 imgfrq -1 ihbfrq 0 ilbfrq 0 iunwrite -1"
            minimize sd 300
            minimize abnr 300
        } until {$ener(GMRS) < 1 }
        my sendCmd "cons harm clear"
        my do {
            minimize sd 300
            minimize abnr 300
            sendCmd "dynamics start verlet nstep 2000 timestp 0.001 nprint 100 -
         iasvel 1 firstt 200 finalt 0 -
         iseed 309234354 -
         inbfrq -1 imgfrq -1 ihbfrq 0 ilbfrq 0 iunwrite -1"
            minimize sd 300
            minimize abnr 300
        } until {$ener(GMRS) < 1}
    }
    if [IOmanipulation::isAAseq $seq] {
        #generate sequence
        #puts "Generate $segid with $seq"
        my sendCmd "read sequ card\n*read from TMD\n*\n [llength $seq] \n $seq \ngenerate $segid setup"
    }
}

oo::define Charmm method build2ndStruct {segid {type alpha} {start 1 } {end end} } {
    set gphi 0;set gpsi 0
    switch  -re $type {
        alpha { set phi -57; set psi -47}
        beta { set phi -125; set psi 122.5}
        loop { set phi [expr rand()*360-180]; set psi [expr rand()*360-180]; set gphi 1;set gpsi 1}
    }
    set nres [my querylength $segid]
    set end [expr "[string map [list end $nres] $end]"]
    for { set i $start} { $i <= [expr $end-2] } {incr i } {
        set j [expr $i +1]
        set k [expr $i +2]
        if $gphi { set phi [expr rand()*360-180] }
        if $gpsi { set psi [expr rand()*360-180] }
        my sendCmd "ic edit\ndihe $i C $j N $j CA $j C $phi\ndihe $j N $j CA $j C $k N $psi\nend"
    }
    my sendCmd "ic param\n"
    if [regexp {.*XMIN = (-?\d+.?\d+)\s*XMAX\s*=\s*(-?\d+.?\d+).*} [my sendCmd "coor stat sele type N end"] matched xmin xmax ] {
        if { [expr abs($xmin -$xmax) ] > 9999.0} { my sendCmd     "ic seed $start N $start CA $start C \n"}
    }
    my sendCmd "ic build\n"
    my sendCmd " mini sd nstep 1000"
}

##Query the length of the segment
oo::define Charmm method querylength { segid } {
    if [regexp {.*>\s+(\d+)\s+atoms have been selected.*} [my sendCmd "coor stat sele type CA .and. segid $segid end"] matched nres] {
        return $nres
    }

}

namespace eval IOmanipulation {
proc isAAseq {seq} {
    set is 1
    foreach aa $seq {
        if ![isAA $aa] {set is 0; break}
    }
    return $is
}
variable AA {ALA AIB ARG ARH ASN ASP CYS ASH ASH1 CYS GLN GLU GLH GLH1 GLY HIS HSD HSC ILE LEU LYS LYH MET PEN PHE PHL PRO SER THR TRP TYR VAL HEME ACE FORM AMN CBX NH2 O2 CO2P CO ETH ST2 TIP3 OH2 SO4 COH ACD ACDH OLA OLAH STE STH SH2 PLM MYR MYRH STHE METM ETP}
proc isAA { aa } {
    variable AA
    set aa [string toupper $aa]
    if { [lsearch $AA $aa]>=0} { return 1}
    return 0
}

}
